#include "../includes/need-a-name.hh"
#include "includes/initcollecter-mhd.hh" 
#include <string>

InitCollecterMHD::InitCollecterMHD(){
  set_parameter("cen1",0.0);
  set_parameter("rhoI",0.0);
  set_parameter("rhoII",0.0);
  set_parameter("mom1I",0.0);
  set_parameter("mom1II",0.0);
  set_parameter("mom2I",0.0);
  set_parameter("mom2II",0.0);
  set_parameter("mom3I",0.0);
  set_parameter("mom3II",0.0);
  set_parameter("b2I",0.0);
  set_parameter("b2II",0.0);
  set_parameter("b3I",0.0);
  set_parameter("b3II",0.0);      
  set_parameter("eI",0.0);
  set_parameter("eII",0.0);
  set_parameter("gamma",5.0/3.0);
  set_parameter("b1zero",0.0);
}

InitCollecterMHD::~InitCollecterMHD(){
}

void InitCollecterMHD::initial_function(double* location,double pillarlen, int currentdim, int index){
  double variable_cache[memory::phynumber];
  double location_cache[gridsetting::dimension];
  double component_cache;  
  double basisInside[memory::dgorder*memory::intergnumber];
  for(int tmpi=0; tmpi<memory::phynumber;tmpi++){
    if(para["sod_shocktube_leve"]==1.0)
      variable_cache[tmpi]=sod_shocktube_leve(location, tmpi);
    if(para["sod_shocktube_stone"]==1.0)
      variable_cache[tmpi]=sod_shocktube_stone(location, tmpi);
    if(para["twoshock_leve"]==1.0)
      variable_cache[tmpi]=twoshock_leve(location, tmpi);
    if(para["briowu"]==1.0)
      variable_cache[tmpi]=briowu(location, tmpi);
    if(para["rj2a"]==1.0)
      variable_cache[tmpi]=rj2a(location, tmpi);
    if(para["rj4d"]==1.0)
      variable_cache[tmpi]=rj4d(location, tmpi);
    if(para["blastwave"]==1.0)
      variable_cache[tmpi]=blastwave(location, tmpi);
  }
  for(int tmpi=0; tmpi<memory::phynumber;tmpi++)
    memory::variable[memory::dgorder*memory::phynumber+tmpi][index]=variable_cache[tmpi];

  if(memory::dgflag==1){
    basis.getBasis(basisInside);
    for(int tmpl=0; tmpl<gridsetting::dimension; tmpl++)
      location_cache[tmpl]=location[tmpl];
    for(int tmpi=0; tmpi<memory::phynumber;tmpi++)
      for(int tmpj=0; tmpj<memory::dgorder;tmpj++){
    	component_cache=0.0;
    	for(int tmpk=0; tmpk<memory::intergnumber; tmpk++){
    	  location_cache[currentdim]=location[currentdim]+(-0.5+tmpk/(memory::intergnumber-1.0))*pillarlen;
	  if(para["sod_shocktube_leve"]==1.0)
	    component_cache+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*sod_shocktube_leve(location_cache,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];	  
	  if(para["sod_shocktube_stone"]==1.0)
	    component_cache+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*sod_shocktube_stone(location_cache,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];
	  if(para["twoshock_leve"]==1.0)
	    component_cache+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*twoshock_leve(location_cache,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];
	  if(para["briowu"]==1.0)
	    component_cache+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*briowu(location_cache,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];
	  if(para["rj2a"]==1.0)
	    component_cache+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*rj2a(location_cache,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];
	  if(para["rj4d"]==1.0)
	    component_cache+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*rj4d(location_cache,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];
	  if(para["blastwave"]==1.0)
	    component_cache+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*blastwave(location_cache,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];
	}
	memory::variable[tmpj*memory::phynumber+tmpi][index]=component_cache;
      }
  }
}

// Figure 14.3 of book: Leveque R. J.
// similar with Fig 8 of Stones J. M. 2008ApJ
double InitCollecterMHD::sod_shocktube_leve(double* location, int phyi){
  double tmp;
  para["cen1"] =0.5;
  if(phyi==0)
    tmp = (location[0]<para["cen1"]) ? 3.0 : 1.0;
  if(phyi==1)
    tmp = 0.0;
  if(phyi==2)
    tmp = 0.0;
  if(phyi==3)
    tmp = 0.0;
  if(phyi==4)
    tmp = (location[0]<para["cen1"]) ? 3.0/(para["gamma"]-1.0) : 1.0/(para["gamma"]-1.0);
  if(phyi==5)
    tmp = 0.0;
  if(phyi==6)
    tmp = 0.0;
  return tmp;
}

// Fig 8 of Stones J. M. 2008ApJ
double InitCollecterMHD::sod_shocktube_stone(double* location, int phyi){
  double tmp;  
  para["cen1"] =0.5;  
  para["gamma"]=1.4;
  if(phyi==0)
    tmp = (location[0]<para["cen1"]) ? 1.0 : 0.125;
  if(phyi==1)
    tmp = 0.0;
  if(phyi==2)
    tmp = 0.0;
  if(phyi==3)
    tmp = 0.0;
  if(phyi==4)
    tmp = (location[0]<para["cen1"]) ? 1.0/(para["gamma"]-1.0) : 0.1/(para["gamma"]-1.0);
  if(phyi==5)
    tmp = 0.0;
  if(phyi==6)
    tmp = 0.0;
  return tmp;  
}

// Figure 14.4 of book: Leveque R. J.
double InitCollecterMHD::twoshock_leve(double* location, int phyi){
  double tmp;  
  para["cen1"] =0.5;  
  para["b1zero"]=0.0;
  if(phyi==0)
    tmp = (location[0]<para["cen1"]) ? 1.0 : 2.0;
  if(phyi==1)
    tmp = (location[0]<para["cen1"]) ? 1.0*3.0 : 2.0*1.0;
  if(phyi==2)
    tmp = 0.0;
  if(phyi==3)
    tmp = 0.0;
  if(phyi==4)
    tmp = (location[0]<para["cen1"]) ? 4.5+1.0/(para["gamma"]-1.0) : 1.0+1.0/(para["gamma"]-1.0);
  if(phyi==5)
    tmp = 0.0;
  if(phyi==6)
    tmp = 0.0;
  return tmp;  
}

// Figure 13 of Stones J. M. 2008ApJ
double InitCollecterMHD::briowu(double* location, int phyi){
  double tmp;  
  para["cen1"] =0.5;  
  para["b1zero"]=0.75;
  para["gamma"]=2.0;
  if(phyi==0)
    tmp = (location[0]<para["cen1"]) ? 1.0 : 0.125;
  if(phyi==1)
    tmp = 0.0;
  if(phyi==2)
    tmp = 0.0;
  if(phyi==3)
    tmp = 0.0;
  if(phyi==4)
    tmp = (location[0]<para["cen1"]) ? 0.5*(1.0+std::pow(para["b1zero"],2))+1.0/(para["gamma"]-1.0) : 0.5*(1.0+std::pow(para["b1zero"],2))+0.1/(para["gamma"]-1.0);
  if(phyi==5)
    tmp = (location[0]<para["cen1"]) ? 1.0 : -1.0;
  if(phyi==6)
    tmp = 0.0;
  return tmp;  
}

// Figure 14 of Stones J. M. 2008ApJ
// RJ2a means fig2a of Ryu Dongsu & Jones TW 1995ApJ
double InitCollecterMHD::rj2a(double* location, int phyi){
  double tmp;  
  para["cen1"] =0.5;    
  para["b1zero"]=2.0/std::sqrt(4.0*constant::dpi);
  para["gamma"]=5.0/3.0;
  if(phyi==0)
    tmp = (location[0]<para["cen1"]) ? 1.08 : 1.0;
  if(phyi==1)
    tmp = (location[0]<para["cen1"]) ? 1.08*1.2 : 0.0;
  if(phyi==2)
    tmp = (location[0]<para["cen1"]) ? 1.08*0.01 : 0.0;;
  if(phyi==3)
    tmp = (location[0]<para["cen1"]) ? 1.08*0.5 : 0.0;;
  if(phyi==5)
    tmp = (location[0]<para["cen1"]) ? 3.6/std::sqrt(4.0*constant::dpi) : 4.0/std::sqrt(4.0*constant::dpi);
  if(phyi==6)
    tmp = (location[0]<para["cen1"]) ? 2.0/std::sqrt(4.0*constant::dpi) : 2.0/std::sqrt(4.0*constant::dpi);
  if(phyi==4)
    tmp = (location[0]<para["cen1"]) ?
      0.5*(std::pow(para["b1zero"],2)+3.6*3.6/(4.0*constant::dpi)+2.0*2.0/(4.0*constant::dpi))+0.5*1.08*(1.2*1.2+0.01*0.01+0.5*0.5)+0.95/(para["gamma"]-1.0) :
      0.5*(std::pow(para["b1zero"],2)+4.0*4.0/(4.0*constant::dpi)+2.0*2.0/(4.0*constant::dpi))+1.0/(para["gamma"]-1.0);
  return tmp;  
}

// RJ4d: Figure 15 of Stones J. M. 2008ApJ
double InitCollecterMHD::rj4d(double* location, int phyi){
  double tmp;  
  para["cen1"] =0.5;  
  para["b1zero"]=0.7;
  para["gamma"]=5.0/3.0;
  if(phyi==0)
    tmp = (location[0]<para["cen1"]) ? 1.0 : 0.3;
  if(phyi==1)
    tmp = (location[0]<para["cen1"]) ? 0.0 : 0.0;
  if(phyi==2)
    tmp = (location[0]<para["cen1"]) ? 0.0 : 0.0;;
  if(phyi==3)
    tmp = (location[0]<para["cen1"]) ? 0.0 : 0.3*1.0;;
  if(phyi==5)
    tmp = (location[0]<para["cen1"]) ? 0.0 : 1.0;
  if(phyi==6)
    tmp = (location[0]<para["cen1"]) ? 0.0 : 0.0;
  if(phyi==4)
    tmp = (location[0]<para["cen1"]) ?
      0.5*(std::pow(para["b1zero"],2))+1.0/(para["gamma"]-1.0) :
      0.5*(std::pow(para["b1zero"],2)+1.0)+0.5*0.3+0.2/(para["gamma"]-1.0);
  return tmp;  
}

// blastwave for 2D/3D MHD test
double InitCollecterMHD::blastwave(double* location, int phyi){
  double tmp;  
  para["rhoI"]=1.0;
  para["pI"]=10.0;
  para["pII"]=1.0;
  para["radius"]=0.15;
  double bzero=1.0/1.415926;
  double distance=0.0;
  double benergy=0.5;
  if(phyi==0)
    tmp = para["rhoI"];
  if(phyi==1)
    tmp =0.0;
  if(phyi==2)
    tmp =0.0;
  if(phyi==3)
    tmp =0.0;
  if(phyi==4){
    for(int tmpi=0; tmpi<gridsetting::dimension; tmpi++)
      distance+=std::pow(location[tmpi]-0.5,2);
    distance=std::sqrt(distance);
    tmp = distance<=para["radius"] ? para["pI"]/(para["gamma"]-1.0)+benergy : para["pII"]/(para["gamma"]-1.0)+benergy;
  }  
  if(phyi==5)
    tmp =bzero;
  if(phyi==6)
    tmp =bzero;
  if(phyi==7)
    tmp =0.0;
  return tmp;  
}
